## ----------------------------- a naive guide for W matrix
## 1. check sparisity: row count and column count of neighbors with positive weights
## 2. check eigen value: the minimum negative eigenvalue should be around -1, 
##                       if the minimum negative eigenvalue is around 0, then the W matrix may not work
## 3. w matrix with complex eigen value (asymmetric w matrix): no worries, just 
##                                                             check 2 for all real eigenvalues



## ----------------------------- output

## [1] "beta"     "rhoy"     "omega"    "sigma2"   "Alpha"   
## [6] "Xi"       "rho0"     "Rho"      "kappa"    "sigma_n2"
## [11] "rhoA"     "L"        "Factor"  

## constant coef
## ar1 parameter for outcome 
## omega for factor loading variance
## error variance
## unit-random coef
## time-random coef
## constant sar coef
## random sar coef
## ar1 parameter for rho_t
## error variance in state space rho_t
## coef for time-vary covar that explains rho_t
## factor loadings matrix 
## factor matrix 

## ---------------------- sim1 block diagnoal w matrix
rm(list = ls())


sink("log/log_block_sim_saom.txt") # start log
begin.time<-Sys.time()


## compile functions 
library(bpNet)
library(mcmcr)
library(coda)
library(MASS)
library(abind)

set.seed(123456789)

## ------------------------- ##
source("code/simulated20529.R")
source("code/net_plot.R")


load("data/SimData/saom_sim1.RData")

mcmc <- 25000
burnin <- 10000

out1.2 <- bpNet(data = data1,   ## a data frame
	              W = W,          ## an array N * N * TT, each slide is a w matrix for a given period
	              index = c("id", "time"), ## id and time
	              Yname = "y",
	              Xname = c("x1", "x2"),
	              Zname = NULL,   ## covar has unit random effect
	              Aname = NULL,   ##           time
	              Contextual = c("x1", "x2"),  ## contectual covar
	              Contextual.effect = "both",  ## contectual effect: "none", "unit", "time", "both"
	              lagY = FALSE,                ## whether to incorporate lagged outcome
	              force = "both",              ## two-way fe" "none", "unit", "time", "both"
	              r = 10,
	              flasso = 1,              ## factor selection
	              rhoZ = as.matrix(rhoZ),  ## time-varying covar that explains rho: T*k     
	              rho_spec = "ar1",                ## random rho: "multilevel", "rw", "ar1"
	              niter = mcmc,
	              burn = burnin)



## Figure A16 plot rho_t
a <- 20
p1.2 <- rho_plot(out1.2, Rho1, constant = FALSE)
p1.2 <- p1.2 +theme(axis.text.x = element_text( hjust = 1, size=a))
p1.2 <- p1.2 +theme(axis.text.y = element_text( hjust = 1, size=a))
ggsave('plots/simsaom/sim_saom_rhot.pdf', p1.2, width = 18, height = 9)


## Figure A15: Factors 
p1 <- omegaPlot(out1.2, oxlim = c(-2, 2), 
	            plotlength = 5, labelname = "Factor")
ggsave('plots/simsaom/sim_saom_factor.pdf', p1, width = 18, height = 20)

print(Sys.time()-begin.time) 
sink() # close log  
